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Summary 

Formal optimization methods and eigenvalue gra- 
dient information are used to develop a stabiliz- 
ing control law for a closed-loop linear system that 
is initially unstable. The method was originally 
formulated by using direct, constrained optimiza- 
tion methods with the constraints being the real 
parts of the eigenvalues. However, because of prob- 
lems in trying to achieve stabilizing control laws, 
the problem was reformulated to be solved differ- 
ently. The method described in this paper uses the 
Davidon-Fletcher-Powell minimization technique to 
solve an indirect, constrained minimization problem 
in which the performance index is the Kreisselmeier- 
Steinhauser function of the real parts of all the 
eigenvalues. The method is applied successfully to 
solve two different problems: the determination of a 
fourth-order control law that stabilizes a single-input 
single-output active flutter suppression system and 
the determination of a second-order control law for 
a multi-input multi-output lateral-directional flight 
control system. Various sets of design variables and 
initial starting points were chosen to show the ro- 
bustness of the method. 

Introduction 

Many optimal control-law synthesis techniques 
require, as a starting point, a control law that results 
in a stable closed-loop system. Throughout this 
paper, such a control law will be referred to as 
a “stabilizing control law.” Sketch A illustrates 
the relationship and sequence of the major steps to 
obtain optimal control laws. 

The methods that are presently used to find 
stabilizing control laws (located in the loop in the 
sketch) are largely trial and error and therefore can 
be time consuming, especially for unstable open-loop 
systems (refs. 1 and 2). An automated method using 
Bass’ algorithm (ref. 3) can determine a stabilizing 
control law when the structure of the control law 
is full-state feedback. However, this method is not 
applicable for reduced-order control laws of arbitrary 
structure. 

The purpose of this paper is to present an au- 
tomated method for obtaining stabilizing reduced- 
order control laws of arbitrary structure. In devel- 
oping this method, an important goal was that these 
control laws be obtained with a minimum amount of 
trial and error by the control-law designer. There- 
fore, the method uses formal optimization methods 
and eigenvalue gradient information to obtain stabi- 
lizing control laws. The method was originally for- 
mulated as a constrained optimization problem with 
the constraints being the real parts of the eigenval- 



Sketch A 


ues. However, because of difficulties in achieving con- 
vergence, the method was reformulated as an indi- 
rect, constrained optimization problem with the per- 
formance index being the Kreisselmeier-Steinhauser 
(KS) function of the real parts of the eigenvalues. 
Analytical expressions for the gradients of the eigen- 
values were derived and implemented. 

The method is applied to solve two example prob- 
lems. The first obtains a fourth-order stabilizing con- 
trol law for a single-input single-output flutter sup- 
pression system. The second determines a second- 
order control law for a multi-input multi-output 
lateral-directional control system. Different choices 
of both the design variables and their initial values 
were used to illustrate the robustness of the method. 

Symbols 

A controller dynamics matrix (M x M) 

B controller input matrix ( M x N 0 ) 

W, Ci type 2 design variables 

C controller output matrix ( N c x M) 

c local chord of wing 

di type 1 design variables (coefficients of 

denominator polynomial) 

F plant dynamics matrix ( N s x N s ) 














F a augmented dynamics matrix (N s + 

M) x (N g + M) 

G u plant input matrix (N s x N c ) 

g % inequality constraints 

H sensor output matrix ( N 0 x N s ) 

I identity matrix 

J performance index 

3 =y/=I 

K matrix equal to real part of [vjU^], 

(N s + M)x {N s + M) 

K 11 (Ng x N s ) matrix formed from 

partitioning matrix K 

K 12 (Ng X M) matrix formed from parti- 

tioning matrix K 

K 21 (M x N g ) matrix formed from parti- 

tioning matrix K 

K 22 (M x M) matrix formed from parti- 

tioning matrix K 

k feedback gain 

fcnom nominal gain 

M order of controller 

N c number of control inputs 

N 0 number of outputs or sensors 

N s number of plant states 

n i type 1 design variables (coefficients of 

numerator polynomial) 

P matrix defined in equation (15a), 

(N c + M)x(N 0 + M) 

Pjk element of matrix P 

R key state selection matrix (M x N s ) 

r drawdown factor 

s Laplace operator 

T(s) transfer function matrix (N c X N 0 ) 

U matrix consisting of all right eigenvec- 

tors, ( N s + M ) x ( N s + M ) 

u control input vector ( N c x 1) 

Uj normalized right eigenvector corre- 

sponding to Aj, (Ng + M) x 1 

V matrix consisting of all left eigenvec- 

tors, (7V S + M) x (N s + M) 


v 2 ; normalized left eigenvector corre- 

sponding to A i, (N s + M) x 1 

X vector of design variables 

Xj jth element of vector X 

x a augmented state vector of order 

(N s + M) x 1 

x c controller state vector of order (Mxl) 

x s plant state vector of order (N a x 1) 

y measurement output vector of order 

(N 0 x 1) 

a iiPiili type 1 design variables 
fyj Kronecker delta 

A* ith eigenvalue 

A i R real part of ith eigenvalue 

Subscripts: 

I imaginary 

R real 

Superscript: 

T transpose 

Matrix notation: 

tr[ ] trace of matrix 

[ ] matrix 

{ } vector 

L J row vector 

[ ] T vector or matrix transpose 

[ ] -1 inverse of square matrix 

Dots over symbols denote differentiation with 
respect to time. 

Equations of Motion in State-Space Form 

The specific applications of the method described 
in this paper are airplane control systems. This 
method is, however, applicable to any linear system 
as long as the equations of motion can be written in 
the following state space form (ref. 2): 

Xg = Fx s + G u u (1) 

y = Hx s (2) 
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These equations are combined with the following 
equations describing a control system with output 
feedback (ref. 2): 


x c = Ax c + By (3) 

u = Cx c (4) 

The block diagram in figure 1 illustrates the con- 
trol scheme that is modeled in equations (1) to (4). 
The dashed lines in figure 1 indicate the three com- 
ponents of the control scheme: the plant, the mea- 
surement, and the control law. Using equations (3) 
and (4), the control law may be expressed in the form 


u = C[sl - A] 1 By (5) 


where the matrix product C[5l — A]” 1 B is the trans- 
fer function matrix T(s). 

Equations (1) to (4) can be combined into a set 
of augmented state equations by first defining an 
augmented state vector as 

<6) 

Next, substituting equation (4) into (1) and equa- 
tion (2) into (3) yields, after rearrangement, 



Defining the augmented system matrix as 


F a = 


F 

BH 


G U C 

A 


(8) 


allows equation (7) to be rewritten as 

x a = F a x a (9) 

Referring back to sketch A in the “Introduc- 
tion,” the initial control law is completely defined 
by matrices A, B, and C in equations (5) and (7). 
Closed-loop stability is determined by evaluating the 
eigenvalues of the augmented system matrix of equa- 
tion (8). The stabilizing control law is completely 
defined by matrices A, B, and C when all the eigen- 
values of equation (8) have negative real parts. The 
method presented in this paper is an automated 
method to select the elements of the A, B, and C 
matrices to provide a stable closed-loop system. 

General Description of Design Algorithm 

The method has been automated and imple- 
mented on the CDC CYBER computer system at the 


Langley Research Center. From now on, this method 
will be referred to as the “design algorithm.” This 
section of the paper presents a general description of 
the design algorithm, shown in flow chart form in fig- 
ure 2. (This flow chart will be referred to frequently 
in this section of the paper.) The purpose of the de- 
sign algorithm is to produce a stabilizing control law, 
indicated by the ellipse on the right side of figure 2. 

The design algorithm employs optimization to 
achieve its purpose. Necessary elements of the op- 
timization routine are the design variables, a per- 
formance index J, and, in the present problem, the 
gradients of the closed-loop eigenvalues with respect 
to the design variables. Beginning at the top of the 
flow chart, the first step in the design algorithm is 
to choose the design variables and their initial val- 
ues, thereby initializing the A, B, and C matrices 
and establishing the initial control law. Next, the 
eigenvalues of the augmented system matrix (closed- 
loop eigenvalues) are computed, after which a deci- 
sion point is reached. If all the eigenvalues have neg- 
ative real parts (that is, are stable), then this initial 
control law is a stabilizing control law and the design 
algorithm is terminated. If one or more eigenvalues 
have nonnegative real parts (that is, are neutrally 
stable or unstable), then the remainder of the design 
algorithm is exercised. The stability of the system 
is checked each time that the eigenvalues are calcu- 
lated. If the system is still unstable, the number of 
iterations is checked to make sure that an arbitrary 
upper limit is not exceeded. This is done to prevent 
the design algorithm from continuing indefinitely if 
the solution is having difficulty converging. 

In the event that the rest of the design algorithm 
will be exercised, a performance index is formulated. 
Gradients of the eigenvalues and the performance in- 
dex are calculated. New values of the design vari- 
ables are chosen by the optimization algorithm and 
the control law is updated. Iteration then continues 
as already described. 

Choice of Design Variables 

This section of the paper deals with choosing de- 
sign variables as indicated by the top box of the flow 
chart in figure 2. The design variables are associated 
with the A, B, and C matrices within the control 
law. The total number of the elements within the 
A, B, and C matrices is M(M + N 0 + N c ). Of this 
total, M(N 0 + N c ) are independent (ref. 4) and can 
therefore be chosen as design variables for the de- 
sign algorithm. The design algorithm accepts design 
variables in either of two different forms. The first 
(referred to a “type 1”) is conveniently understood 
and interpreted by use of a transfer function matrix 
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and concepts from classical control theory; the sec- 
ond (type 2) is similar to a linear quadratic Gaussian 
(LQG) formulation. Both types of design variables 
result in reduced-order control laws. 

Type 1 design variables. Type 1 design variables 
are chosen to be specific elements of the A, B, and 
C matrices. Any specific transfer function can be 
realized into state space A, B, C form in a variety 
of ways. For the sake of illustration, let us assume 
a fourth-order control law (M — 4) for a single- 
input ( N c = 1) single-output (N a = 1) system. The 
number of available design variables for the design 
algorithm is M(N a + N c ) - 8, and the transfer 
function matrix reduces to a single transfer function 
that may be written as 


T(a) = C[sl - A] - *B = . ^ 3 « 3 + » 2 a 2 +ni a + n 0 ) 
-f- + d\s + dQ 


(10) 


A state space representation of the control law shown 
in equation (10) is given as follows: 



-0 0 0 

-d o' 

A = 

1 0 0 

~di 


0 1 0 

-d 2 


.0 0 1 

-do J 


r«°i 


B = < 

i 1 \ 



n 2 ( 



1 no ) 


o 

o 

o 

II 

O 



(11a) 

(lib) 

(11c) 


Within the right side of equation (10), the nine quan- 
tities (fc(Gain); no, ni, n 2 , and n 3 (coefficients of the 
numerator polynomial); and do, d\, d 2 , and do (co- 
efficients of the denominator polynomial)) are read- 
ily identified as elements of the A, B, and C matri- 
ces and are therefore candidates for design variables. 
From these nine quantities, the control-law designer 
chooses up to eight. 


Type 2 design variables. Assuming again a 
fourth-order control law for a single-input single- 
output system, there are again eight available design 
variables. The type 2 design variables are the four 
elements of the B matrix and the four elements of 
the C matrix where 


f M 



Kht 

C = (.CO Cl c 2 C 3 J 


(12a) 

(12b) 


These design variables are analogous to elements of 
the gain matrices obtained from the LQG solution. 

Recalling that the control laws in this paper are 
of reduced order, a key state selection matrix R 
is employed for retaining, deleting, and combining 
plant states. With the B and C matrices given in 
equations (12a) and (12b), respectively, the A matrix 
is a function of the B and C matrices and is given 
by (ref. 2) 


A = RFR t - BHR r + RG W C (13) 


Eigenvalues of Augmented System Matrix 

Once the control law has been initialized, the aug- 
mented system matrix is calculated. The eigenvalues 
of this matrix are calculated by subroutine EIGEN 
in ORACLS (ref. 5). The stability of the system is 
checked; if all the real parts of the eigenvalues are 
negative, then the program is terminated. If not, an- 
other decision must be made. If a maximum limit on 
the number of iterations is reached the algorithm is 
terminated. If not, the rest of the design algori thm 
is executed. 

Eigenvalue Gradients 

The eigenvalue gradients are needed for the de- 
sign algorithm. Analytical expressions for the eigen- 
value gradients are derived and depend on the choice 
of the design variables. Although the complete 
derivation of the expressions is given in appendix A, 
a general definition of the gradients of the eigenvalues 
is given by 


$ = [G' T v i ufH' T ] 


where 


P = 


0 C 

B A 


\{N c +M)x{N 0 +M) 

Matrices G / and H / are defined, respectively, 


G' = 

H' = 


Gr« 0 

0 I 

H 0 

0 I 


(iy s +M)x(AT c +M) 

(N 0 +M)x(N a +M) 


(14) 


(15a) 

by 

(15b) 

(15c) 


where I is an ( M x M) identity matrix. 

As seen in equation (14), both the left and right 
eigenvectors (v* and u;, respectively) are needed 
to calculate the eigenvalue gradients. The right 
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eigenvectors were calculated by solving the following 
eigenvalue problem: 

Auj = (16) 

There are two ways of calculating the left eigen- 
vector. The first method is to solve the eigenvalue 
problem shown in equation (16) again but substi- 
tuting A t in place of A. The second method is to 
calculate the left eigenvector directly from the right 
eigenvector. This method is derived from the relation 


type 2 design variables are 

= [K 2 i]H t - [K 22 ]RH r (20a) 

= G% [K 12 ] + G£R T [K 22 ] (20b) 

Optimization Formulation 


vjuj = 6ij (17) 

where 6{j is the Kronecker delta that equals 1 when 
i = j and equals 0 when i ^ j. For a given eigenvalue 
with its corresponding normalized eigenvectors, the 
inner product of the left and right eigenvectors is 
unity. For normalized eigenvectors corresponding to 
different eigenvalues, this inner product is 0. 

The second method mentioned previously can be 
implemented in different ways. One way is to invert 
the complex matrix formed by all the right eigen- 
vectors. To reduce the computation time, another 
method was developed. In this method, the left 
eigenvectors were calculated directly from the right 
eigenvectors without using complex matrix opera- 
tions. A derivation of the method is given in ap- 
pendix B. 

The expressions for the eigenvalue gradients de- 
pend on the choice of the design variables. For type 1 
design variables, the expressions for the gradients are 
given by 

= [K 22 ] (18a) 

= [K 2 i]H t (18b) 

^ = g£[K 12 ] (18c) 

where [K22 ]i [K21], and [K12] are partitioned matri- 
ces from the real part of v^uj 1 given as 


L L X *-21J L iV 22J J (jV s +M)x(JV a +M) 

(19) 

For type 2 design variables, the additional relation- 
ship given by equation (13) must be satisfied. Be- 
cause of this interdependence, the gradients of the 
eigenvalues with respect to the type 2 design vari- 
ables have an additional term. The expressions for 
the gradients of the eigenvalues with respect to the 


Direct constrained approach . A stabilizing con- 
trol law requires that the real parts of all closed-loop 
eigenvalues be negative. During the early develop- 
ment of the design algorithm, it seemed appropri- 
ate to meet this requirement directly by the use of 
constrained optimization. Therefore, inequality con- 
straints of the form 

&<0 (* = 1 , .... n) (21) 

were specified where are the real parts of the 
eigenvalues. 

For any optimization problem, an objective func- 
tion to be minimized must be specified. Notice that 
the purpose of the design algorithm will have been 
met when all the constraints are satisfied. There- 
fore, the selection of the objective function for this 
formulation is only of secondary importance and it 
was arbitrarily chosen to be one-half the sum of the 
squares of the elements of the row vector C: 

J = l -CC T (22) 

The CONMIN program (ref. 6), employing the 
method of feasible directions (ref. 7), was chosen for 
incorporation in the design algorithm. 

An insurmountable numerical problem occurred 
within the design algorithm that resulted in the aban- 
donment of the direct optimization approach. The 
subroutine that calculates the eigenvalues sorts them 
in order of increasing magnitude. During the opti- 
mization process, the eigenvalues do not remain in 
the same order but they may switch places. Unsta- 
ble eigenvalues may switch places with stable ones, 
thus making it very difficult to keep track of unstable 
eigenvalues. Because gradients of the violated con- 
straints (unstable eigenvalues) are needed during the 
optimization process, and since it is difficult to track 
the eigenvalues, it appears to the optimization rou- 
tine that the gradients of the unstable eigenvalues are 
discontinuous. Therefore, as stated previously, this 
formulation was abandoned. 
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Indirect constrained approach. The problem was 
reformulated by using indirect constrained optimiza- 
tion methods. The performance index was chosen 
to be the Kreisselmeier-Steinhauser (KS) function 
(ref. 8). The KS function, which is a penalty func- 
tion that takes into account the real parts of all of 
the eigenvalues (stable and unstable), is given as 



In 


n 

Y2 ex p( rA * ft ) 


u=l 


(23) 


The function is dominated by the unstable eigenval- 
ues (those with positive real parts). The contribu- 
tion of the stable eigenvalues to the KS function is 
relatively small because of the behavior of the ex- 
ponential function. The quantity r is a drawdown 
factor that (for values of r greater than unity) serves 
to increase the relative contribution of the unstable 
eigenvalues to the performance index. Figure 3 shows 
the exponential function and its relation to some rep- 
resentative eigenvalues. 

The gradients of the performance index with re- 
spect to the design variables were derived analyti- 
cally. The equations for the gradients are 


<9KS 

dX 3 - 


E exp{r\ iR )dX iR /dX j 
i = 1 

n 

E ex P(^; fi ) 

i = 1 


(24) 


Expressions for the partial derivatives dX iJdXj 
within the numerator of equation (24) are given in 
equations (18) and (20). 

To perform the indirect optimization, the Auto- 
mated Design Synthesis (ADS) package of programs 
(ref. 9) was used. Within ADS, a variable met- 
ric method, the Davidon-Fletcher-Powell algorithm 
(ref. 10), was chosen because of its known rapid rate 
of convergence for this type of problem. 

Termination 

The design algorithm may terminate for the fol- 
lowing three different reasons, as indicated in fig- 
ure 2. 


Successful termination. The first reason is the 
successful termination that occurs when a stabilizing 
control law is found. 

Maximum iterations exceeded. The second rea- 
son occurs when the maximum number of iterations 
is exceeded. Depending on how the optimization pro- 
cess is proceeding, a feasible solution may be reached 


if one just increases the maximum number of itera- 
tions or restarts the design algorithm at the final con- 
trol law of the previous optimization process. If, in 
the attempt to minimize the performance index, the 
magnitude of the positive real part of the unstable 
eigenvalues continues to decrease, but the number of 
unstable eigenvalues starts to increase, then a solu- 
tion may not be achievable. In this case, the only 
option is to choose another initial starting point. 

Unsuccessful termination. The third reason oc- 
curs when a local minimum of the performance index 
is found before the closed-loop system is stabilized. 
In this case it was postulated that convergence could 
be reached if the design algorithm started with this 
control law for the next iteration but with an in- 
creased value of the drawdown factor r. An option 
to restart the optimization process with an increased 
value of r was implemented in the design algorithm; 
but after exercising the algorithm with several dif- 
ferent starting points, it was determined that if the 
algorithm did not converge with the initial value of 
r, it would not converge with an increased value of 
r. Therefore, the option was eliminated. The default 
value of r used in all cases was chosen as unity. How- 
ever, this value is checked initially to make sure that 
the magnitude of r multiplied by the value of the real 
part of the most unstable eigenvalue did not exceed 
a maximum value. This maximum value is defined 
by the computer and represents an upper limit for 
which the exponential function can be evaluated. In 
a few test cases run from different starting points, a 
local minimum was found before a feasible solution 
was found. 

Numerical Examples 

The design algorithm described in this paper is 
applied to solve two different problems. The first ap- 
plication determines a fourth-order control law that 
stabilizes a single-input single-output unstable open- 
loop system for active flutter suppression. The sec- 
ond application determines a second-order control 
law for a multi-input multi-output, open-loop unsta- 
ble lateral-directional control system. Within each 
application several examples will be given to illus- 
trate the effectiveness of the method. 

First Application 

The configuration used for the flutter suppression 
problem is an aeroelastic wind-tunnel model of the 
DAST ARW-1 (Drones for Aerodynamic and Struc- 
tural Testing— Aeroelastic Research Wing 1) shown 
in figure 4. The plant model is a single-input single- 
output system of 25th order. The locations of the 
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sensor and the control surface are shown in the fig- 
ure. The sensor measures the vertical acceleration 
that is fed back through the control system, which 
in turn determines control deflections needed to sup- 
press flutter. 

Case 1, The first example involves only one de- 
sign variable: the gain of a given fourth-order control 
law. An unstable starting point was derived from a 
known stabilizing control law, given by equation (30) 
in reference 2, and is repeated here as 

k{s + 24.74)(s 2 + 87.63s + 13 8 °6) 

^ “ (s + 3.864) (s + 3270) (s 2 + 20.97 s + 1423) 

(25) 

When the gain k is the nominal value of gain in 
reference 2 (fc n om)» this control law results in a stable 
closed-loop system. A value of gain (k = 10fc n om) 
was chosen to provide an unstable closed-loop system 
as a starting point. Using the design algorithm, a 
stable closed-loop system was found. 

The history of the normalized performance index 
is given in figure 5. The solution converged to 
a stabilizing control law during the fifth iteration. 
The final value of the normalized performance index 
was —0.044. The actual value is not important 
and can be positive, negative, or zero. For this 
particular example with only one design variable, the 
optimization algorithm corresponded to a simple one- 
dimensional search. Figure 6 shows a history of the 
gain during the iteration process. The final value of 
gain is 2.08 k nom . Figure 7 shows a history of some 
of the eigenvalues. At the starting point, there was 
one complex conjugate pair of unstable eigenvalues 
indicated by the open-circle symbol in the right half- 
plane. 

Case 2. The second example uses the same form 
of the control law as in equations (10) and (25): 
a third-order numerator polynomial over a fourth- 
order denominator polynomial. The design variables 
were chosen as the coefficients of the numerator and 
denominator polynomials (type 1 design variables). 
The initial values of the eight design variables as well 
as the value of gain were all chosen as unity. 

Figure 8 shows a history of the normalized per- 
formance index. The design algorithm converged to 
a stable system in five iterations. Figure 9 shows 
histories of the design variables. Design variables 
do, di, and no (corresponding, respectively, to the 
coefficients of the two lowest order terms of the de- 
nominator polynomial and the constant term of the 
numerator polynomial) remained at their initial val- 
ues. A history of some of the eigenvalues is given in 


figure 10. At the starting point, there were a com- 
plex conjugate pair of unstable eigenvalues and a real 
unstable eigenvalue. During the first iteration the 
real eigenvalue became more unstable. During the 
second iteration it combined with a real stable eigen- 
value and formed a complex conjugate pair of unsta- 
ble eigenvalues. The design algorithm then stabilized 
this and the other unstable pair of eigenvalues. 

Case 3 . The method was also verified by using 
type 2 design variables. The order of the control 
law was chosen again to be four, resulting in eight 
design variables (four elements each of the B and C 
matrices). To generate the A matrix of the control 
law, four states needed to be chosen. The key states 
chosen were the displacements and velocities of the 
first wing bending and first wing torsion modes. 

An arbitrary starting point was chosen in which 
all the design variables were set to unity. However, 
for this starting point, the design algorithm did not 
converge to a stabilizing control law because a local 
minimum was found without obtaining a stabilizing 
control law. A new starting point was arbitrarily 
chosen by setting the elements of the C matrix to 
0.1 (an effect which is equivalent to decreasing the 
gain of the control law by a factor of 10). From this 
starting point, the design algorithm converged to a 
stabilizing control law in 14 iterations. A history 
of the normalized performance index is shown in 
figure 11. Histories of the design variables are shown 
in figure 12. All the values of the design variables 
changed during the optimization process. A history 
of some of the eigenvalues is shown in figure 13. 
At the starting point, there were a pair of complex 
conjugate unstable eigenvalues and two real unstable 
eigenvalues. 

Three observations about the behavior of the 
eigenvalues can be made from figure 13. The first 
is that during the attempt to stabilize all the eigen- 
values, the eigenvalues can actually move in the di- 
rection opposite to that desired: that is, unstable 
eigenvalues become more unstable and stable eigen- 
values become less stable. The second is that during 
the optimization process, some eigenvalues can criss- 
cross back and forth across the imaginary axes, thus 
becoming unstable or stable until all the eigenvalues 
are finally stable. The third is that real eigenval- 
ues combined and formed complex conjugate pairs of 
eigenvalues. 

Second Application 

The second application successfully solved by us- 
ing this method was the determination of a second- 
order control law for stabilizing an unstable multi- 
input multi-output lateral-directional automatic 
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flight control system (AFCS) for the DAST ARW -2 
(ref 11 ). The augmented ARW -2 vehicle has an un- 
stable Dutch roll mode. Figure 14 shows a planform 
of the ARW -2 with the locations of the control sur- 
faces (rudder and differential stabilizer) and the sen- 
sors (which measure roll and yaw rates) to be used in 
the lateral-directional AFCS. Figure 15 is a block di- 
agram of the lateral-directional AFCS, which is char- 
acterized by two feedback loops (one for yaw rate and 
one for roll rate) and cross-feed. The design algo- 
rithm was used to find the coefficients of the two first- 
order transfer functions identified as control laws in 
the block diagram. These transfer functions from 
figure 15 are given as 


t *« 

(26) 


(27) 


The transfer function is denoted by Ti(s) in the roll 
rate loop and by T 2 (s) in the yaw rate loop. The 
transfer function matrix representation of the control 
law can be realized in the state space A, B, C form 
given by the 2 x 2 matrices shown in 


A = 
B = 

C = 


— <*1 0 

0 —02 

'01 o' 

0 02 . 

'll o ' 

. 0 72. 


(28) 

(29) 

(30) 


Two different examples are presented, both with six 
type 1 design variables. The six design variables are 
a i, a 2, 01, 02, ll, and 72 . 


Case 1. For the first example, the initial values of 
the design variables in the order listed previously are 
0, 1, 1, 1, 1 , and 1, which correspond to Ti(s) = 1 /s 
and T 2 (s) = l/(s + 1 ). A history of the normal- 
ized performance index is shown in figure 16. The 
program converged to a stable system in one iter- 
ation. The design variables are shown in figure 17. 
All the design variables changed during the optimiza- 
tion. A history of some of the eigenvalues is shown 
in figure 18. At the starting point there was a pair 
of unstable complex conjugate eigenvalues that were 
stabilized by using the design algorithm. 


Case 2. For the second example, the initial values 
of the design variables are 0 , 0 , 1 , 1 , 0 , and 0 , which 
correspond to T^s) = 0/s and T 2 (s) = 0/s. The 
eigenvalues at this starting point are the open-loop 
eigenvalues plus two additional eigenvalues at the 
origin. A plot of the normalized performance index 
is shown in figure 19. A stable system was found in 
three iterations. The histories of the design variables 
are shown in figure 20. A history of some of the 
eigenvalues is shown in figure 21 . 

Concluding Remarks 

A method has been developed that uses a formal 
optimization method and eigenvalue gradient infor- 
mation to determine a control law that stabilizes a 
linear system. The goal of the method is to determine 
a stabilizing control law of a given general structure 
with a minimum amount of trial and error by the 
control-law designer. The method is easy to imple- 
ment for a general problem in which the state space 
model of the system and the structure of the control 
law are given. The method was originally formulated 
as a constrained optimization problem, but because 
of convergence problems, the problem was reformu- 
lated as an indirect, constrained optimization prob- 
lem in which the performance index is a function of 
the real parts of the eigenvalues. Analytical expres- 
sions for the gradients of the eigenvalues were derived 
and implemented. 

The effectiveness of this method was successfully 
demonstrated by using several numerical examples. 
Using a single-input single-output flutter suppression 
problem as the first application of the method, sev- 
eral fourth-order control laws were determined that 
stabilized the unstable system. Two types of the de- 
sign variables and several initial values of these de- 
sign variables were used. The second application was 
the determination of a second-order control law for a 
multi-input multi-output lateral-directional system. 
The results of the numerical examples showed that 
the method was very successful in reaching its objec- 
tive of finding a stabilizing control law. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
August 13, 1985 
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Appendix A 
Eigenvalue Gradients 

This appendix describes the derivation of the 
equations for the eigenvalue gradients (shown in 
eqs. (18) and (20)) using the method of reference 12. 
It is known that for any real square matrix A, that 
is ( N x IV), the following equations are satisfied: 

Aui = (Al) 

A T Vj = A {Vi (A2) 

In the previous equations, A; are distinct eigenvalues, 
and u i and v t are complex right and left normalized 
eigenvectors (dimensioned N x 1) associated with 
the specific eigenvalue. Differentiating equation (Al) 
with respect to a parameter p gives 


dA 


du,- d A. 


duj 


— — u, : + A-r-^ = -^rrUj -t- A*-^— (A3) 


dp 


dp dp 


dp 


Equation (A3) is a vector equation. In order to 
solve for dA Jdp, which is a scalar, the equation is 
rearranged and premultiplied by to obtain 


'dA*. 

dp 


T dA 

- v * dT* 


+ V. 


T A dUj 

dp 


-yj\ 


dui 

dp 


(A4) 


Since this is a scalar equation, the last two terms can 
be transposed since the transpose of a scalar is itself. 
Using the distributive property 


dp 


v r^i u .- v T^ U ' + (7Ui T 

v » — u * - v * g p u * + 


— ( A T Vj - AjVj) (A5) 


dp 


The last term is identically 0 by employing equa- 
tion (A2). Equation (A5) becomes 


’dA i 

dp 


-dA 


Ui v * dp Ui 


(A6) 


The parameters p are defined to be the elements Pj ^ 
of the matrix P where P is defined as 


P = 


0 C 

B A 


(Nc+M)x{N 0 +M) 


(A9) 


N c is the number of control inputs, N 0 is the number 
of outputs, and M is the order of the controller. 

The matrices F', G', and H' are defined, respec- 
tively, as 


-l-l/ 

F 

0 

F = 

.0 

0 

G' = 

G u 

0 

( 

] 

H' = 

H 

0 

0 

I 


(N a +M) x (N a +M) 

r 

. ( N a +M)x(N c +M ) 
{N 0 +M)X(N S +M) 


(A10) 

(All) 

(A12) 


The symbol I is an (M x M) identity matrix. Using 
the definitions given in equations (A10), (All), and 
(A12), the augmented state system matrix is then 
defined by 

F a = F' + G'PH' (A13) 

Differentiating equation (A13) and noting that F' is 
not a function of the parameters Pjk, the following 
equation is obtained: 


dF a d(G'PH') 
dPjk ~ dP jk 


(A14) 


Substituting equation (A14) into (A7) gives the fol- 
lowing scalar equation: 

(A15) 

dPjk * dPjk 

The trace properties of compatible matrix products 
can be used to write equation (A15) as 


Since dAj/dp is a scalar, equation (A6) is a scalar 
equation, and vfu* = 1 for normalized eigenvalues, 
the following equation is obtained: 


dXi 

dp 



(A7) 


Steps (Al) to (A7) are also given in reference 12. 

The next effort is to evaluate dA/dp. Matrix A 
is the closed-loop augmented system matrix. For our 
specific example, the augmented system matrix A is 
the matrix F a defined by 


F a 


F 

BH 


G U C 

A 


(N s +M)x(N 3 +M) 


(A8) 


dXj 

dPjk 


= tr 


d(G'PH') 

dPjk 


u *v- 


(A16) 


Since H' is not a function of Pj k , equation (A16) can 
be written as 


dXj 

dPjk 


= tr 


d(G'P) 

dPjk 



(A17) 


Since tr[AB] = tr[BA], equation (A17) can be 
written as 


dXj 

dPjk 


= tr 


T^(G'P) 

HUiVi ~dP^r 


(A18) 
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Since tr[A] — tr[A r ], equation (A18) can be written 
as 


Substituting equation (A25) into (A21) gives 


dXi 


dP, 


= tr 


jk 


a ( pTG ' T ) T„,T 


or 


jk 


y i u i H 


(A19) 


Since G' is not a function of Pj k , equation (A19) can 
be written as 


d\i 


BP, 


= tr 


jk 


d{P T ) fT 


dP„ 


G ViU 


jk 


t ufH ,T = f 


G ,T VjufH ,T 


jk 

(A20) 


Combining all derivatives with respect to Pj k gives 


§ = [G' r v(ufH 


*71 


(A21) 


Having already defined G' and H' in equa- 
tions (All) and (A12), respectively, the next step 
is to get an expression for v*uf . The right eigenvec- 
tor Uj is calculated by solving the eigenvalue problem 
given in equation (Al). The left eigenvector v 4 - is cal- 
culated from the right eigenvector without using the 
complex matrix operations. (This method is derived 
in appendix B.) The product of the left and right 
eigenvectors is expanded by separating both the left 
eigenvector and the right eigenvector into their re- 
spective real and imaginary parts as seen in 


y i u i = { y i R + 3 y ii) ( u i R + J u i,) T (A22) 
Expanding the product of the eigenvectors gives 
Vi uf = (v iR uf R - v^ uf,) +j (y iR ul + v t ,uf fl ) 

(A23) 

For the gradients of the real parts of the eigenvalues, 
only the real part of equation (A23) is needed. This 
term is shown as 


K = Real [v*ufj (A24) 

where K can be partitioned as shown in 


[K] = 


'Ku k 12 1 


cs 

<N 

<N 

1 



(N a x N s ) 
(M x N a ) 


{N s x M) 
(M x M) 


(A25) 


dX i 

dP 


Git 0 

T 

Ku 

Ki 2 

'h o' 

0 I 


k 21 

k 22 

0 I 


T T 


(A26) 


Knowing the definition of P given in equation (A9), 
the gradients of the real parts of the eigenvalues 
with respect to the matrices A, B, and C are given, 
respectively, by 


a? = PM (A27) 

^ = [K 21 ]H t (A28) 

^ = G^[K 12 ] (A29) 

These equations are valid for type 1 design variables. 

For the case 2 design variables, the following 
equation is also satisfied: 

A = RFR t - BHR r + RG U C (A30) 

The total derivatives of the eigenvalues with respect 
to the design variables B and C have an additional 
term because of this expression. The gradients of 
the real parts of the eigenvalues for type 2 design 
variables are shown as 


n H r 

dB dB dA 

, cT-pT ^jiR 

dC dC ^ M dA 


(A31) 

(A32) 


After substituting the appropriate partial derivatives 
(eqs. (A28) and (A29)) into equations (A31) and 
(A32), the gradients of the real parts of the eigen- 
values with respect to the type 2 design variables are 
given as 


= [K 2 i]H t - [K 22 ]RH t (A33) 

= G£ [K 12 ] + G„R r [K 22 ] (A34) 
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Appendix B 


Since the imaginary part of the right eigenvector is 
0, the following equations are satisfied: 


Calculation of Left Eigenvector 

This appendix describes how the complex left 
eigenvectors are calculated directly from the com- 
plex right eigenvectors without complex matrix op- 
erations. The definitions of the right eigenvectors u, t 
and left eigenvectors v* are given, respectively, by 

Auj = AjUj (Bl) 

A T Vj = AjVj (B2) 


y ii u 3i = 0 

v L u n = 0 


(BIO) 

(BH) 


For a right eigenvector uy that is complex, the 
defining relations are equations (B6) and (B7). Its 
complex conjugate is also an eigenvector. Equa- 
tion (B7) becomes 


-yT u . -)- y T u . — 0 

*r3fi 


(B12) 


Note that both A and A T have the same eigenvalues 
and if A = A r , the right and left eigenvectors are 
equal. The following equation is known to be satisfied 
for normalized eigenvectors: 

yfuj = Sij (B3) 

Since the eigenvectors u j and Vj are both complex, 
equation (B3) can be expanded to 

(v<* + JVij) { u m + 3 u h) = 6 ij ( B4 ) 

Expanding the relation, the equation becomes 

( y L u m- y ii u h) + j ( v i R u h +y ii u m) = 

(B5) 

The real terms can be separated from the imaginary 
terms to form 


y L U jR- y l U 3l= 8 ^ 

y iR U 3l +y il U jR=° 


(B6) 

(B7) 


If the right eigenvector is real where Uj is given 
and i refers to any eigenvector, equations (B6) and 
(B7) become, respectively, 


Adding equations (B12) and (B7) gives 

2 ( v 5 u m)=° 

which can also be written as 


v T u . — — \T w . — o 

V n U JR V *J JR 


(B13) 


(B14) 


Equation (B14) can be substituted into equation (B7). 
Then, equations (B6) and (B7) can be added to 
obtain 


v iR U m - v 5 u i/ + y iR U h - y l n 3R = 6 v (B15) 
Applying the distributive property gives 


( V *R ~ y il) T ( U 3R + U 3l) = S i3 ( B16 ) 

Combining all the equations for all the eigenvalues 
gives 

(Vi ? -V/) T (U fl + U / ) = I (B17) 

The matrix of the left eigenvectors can be determined 
from the right eigenvectors by using 


(Vfi-V/)= [(U^ + U/)- 1 


(B18) 


y L U 3R = b i3 
y il U 3R=° 


(B8) 

(B9) 


The real and imaginary parts of the eigenvectors 
corresponding to each specific eigenvalue can then 
be extracted. 
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Figure 2. Flow chart of design algorithm. 
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Figure 8. History of normalized performance index for first application — case 2. Normalization factor, 339.5172. 
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Figure 11. History of normalized performance index for first application — case 3. Normalization factor, 66.954. 
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Figure 13. History of eigenvalues for first application — case 3. 










Figure 14. DAST ARW-2 sensor and control-surface locations for lateral-directional automatic flight control 
system (AFCS). 
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Figure 17. History of design variables for second application — case 1. 








0 12 3 


Number of iterations 

Figure 19. History of normalized performance index for second application — case 2. Normalization factor, 
1.6898. 
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